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Abstract 

We focus on two approaches that have been proposed in recent years for the explanation of the 
so-called FPU paradox, i.e. the persistence of energy locahzation in the 'low-g' Fourier modes of 
Fermi-Past a~Ulam nonlinear lattices, preventing equipartition among all modes at low energies. 
In the first approach, a low-frequency fraction of the spectrum is initially excited leading to the 
formation of 'natural packets' exhibiting exponential stability, while in the second, emphasis is 
placed on the existence of '(/-breathers', i.e periodic continuations of the linear modes of the lattice, 
which are exponentially localized in Fourier space. Following ideas of the latter, we introduce in this 
paper the concept of 'g-tori' representing exponentially localized solutions on low-dimensional tori 
and use their stability properties to reconcile these two approaches and provide a more complete 
explanation of the FPU paradox. 
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I. INTRODUCTION 

In a number of recent papers, Flach and co-workers [sl, Q] discussed the role of simple 
periodic solutions, called 'g-breathers', in the dynamics of the a— and /5— versions of the 
Fermi-Pasta-Ulam model with fixed boundaries. A g-breather is the continuation, for a 7^ 
or /? 7^ 0, of the simple harmonic motion exhibited by the g— th-mode in the uncoupled case 
(a = /3 = 0) and implies motion by a unique frequency, which is nearly equal to the frequency 
of the g-th mode in the uncoupled case. For small values of the coupling parameters (and 
the energy), the distribution of energy in g-breathers stays localized in practice upon only 
a few modes. Thus, g-breathers offer new insight in understanding the problem of energy 
localization, as well as the long term deviations from equipartition of the energy among the 
modes, i.e. the origin of the FPU paradox ^. (For a recent and comprehensive review on 
discrete breathers see js]). 

Analytical estimates of various scaling laws concerning g-breathers can be obtained via 
the method of Poincare-Lindstedt series jll, Q] ■ In particular, expanding the solutions for all 
the canonical variables up to the lowest non-trivial order with respect to the small parameters 
yields an exponentially decaying function for the average harmonic energy of the gth mode, 
i?(g) oc exp(— 6g). The value of h depends on: i) a and/or /5, ii) the number of degrees of 
freedom A^, and iii) the total energy E given to the system. The most important property 
of the g-breathers is that their energy profiles -E'(g) are quite similar to those of 'FPU- 
trajectories', i.e. solutions generated by initially exciting one or a few low-g modes, whose 
exponential localization has been noted since a long time ago , lOj. On the basis of 

this similarity it has been conjectured that there is a close connection between g-breathers 
and the energy localization properties of FPU-trajectories. It is intriguing, however, that 
this localization persists even for values of the parameters (coupling, energy and A^) for 
which the corresponding g-breather solution has become unstable. A heuristic argument for 
interpreting this phenomenon was offered by Flach and Ponno 

Adopting a different approach, Berchialla et al. [l^ explored in detail the question of 
energy equipartition in FPU experiments where a constant (low-frequency) fraction of the 
spectrum (instead of just one mode) was initially excited. The numerical indication from this 
study, based on a limited range of values of E and A^, was that the flow of energy across the 
high-frequency parts of the spectrum takes place exponentially slowly, by a law of the form 
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T oc exp(e~^/^), where T is the time needed for the energy to be nearly equally partitioned, 
and 6 = E/N is the specific energy of the system. In fact, this dependence appears as a 
piecewise power law, with different 'best-fit' slopes in different ranges of values of e. For 
example, a power-law behavior of the type T oc e^^ was found in a subinterval of e values 
considered in [l^, which fits nicely previous results on the scaling of the equipartition time 
with e beyond a critical 'weak chaos' threshold reported in ll||. Nevertheless, the question 
of whether the scaling laws characterizing the approach to equipartition depend on the 
total energy E, or the specific energy E/N is still open, since no rigorous results have been 
provided so far in the literature. On the other hand, numerical results are available over 
a limited range of values of E and A^, in which E is varied proportionally to N. Various 
semi-analytical or numerical approaches to this question are reviewed in 12|. 

Interestingly enough, even if one starts by exciting one mode with large enough energy, 
one observes, long before equipartition, the formation of metastahle states coined 'natural 
packets' 131], in which the energy undergoes first a kind of equipartition among a group of 
low-frequency modes, as if only a fraction of the spectrum was initially excited. We may 
thus conclude that the phenomenon of metastability characterizes FPU-trajectories resulting 
from all types of initial excitations of the low frequency part of the spectrum (see for 
a review of the 'history' of the metastability scenario in the FPU problem). Furthermore, 
empirical scaling laws can be established |13l] concerning the dependence of the 'width' of 
a packet on the specific energy e. These are consistent with the laws of energy localization 
obtained via either a continuous Hamiltonian model which interpolates the FPU dynamics 



in the space of Fourier modes 
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16|, or the g-breather model A difference, however. 



between the two models is that in the framework of the continuous model the constant b in 
E{q) oc exp(— 6g) turns out to be independent of E. 

The results reported in the present paper aim to provide a more complete explanation 
of the FPU paradox of energy non-equipartition by reconciling the presence of g-breathers 
and their induced energy localization on the one hand, with the occurrence of metastable 
packets of low-frequency modes on the other. 

To extend the results obtained for FPU-trajectories in 10, lUl], let us observe that the 
packet of modes excited in these experiments corresponds to the modes with spectral num- 
bers satisfying the condition (A^ + l)/64 < g < 5(A^ + l)/64, with N equal to a power of 
two minus one. For simplicity, let us alter this slightly and consider, instead, the condition 
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1 < g < 4A^/64, with a power of two. The lowest possible value of allowed is = 16, 
for which the above condition implies that only the q = 1 mode is initially excited, giving 
a solution close to a g-breather, which is an orbit lying on an invariant one-torus of the 
system. Now, if is doubled (A^ = 32), the same condition implies that modes q = 1 and 
q = 2 may now be excited, meaning that the resulting FPU-trajectories may be regarded 
as lying close to an invariant two-torus of the system A^ = 32. In general, for A^ = 16s, the 
modes g = 1, 2, . . . , s are initially excited and the motion should be regarded as lying close 
to an invariant s-torus of the system. The same holds true when natural packets of width s 
are formed in experiments in which the initial conditions are as adopted in 13l |. 

This leads to the idea that the properties of the FPU-trajectories could be understood 
by considering classes of special solutions lying not only on one-dimensional tori (as is 
the case with g-breathers) , but also on tori of any low dimension s << N, i.e. solutions 
with s independent frequencies, representing the continuation of motions resulting from 
exciting s modes of the uncoupled case. Generalizing the concept of g-breathers, we call 
q-tori the quasiperiodic solutions solutions on such low-dimensional tori. The main body 
of the present paper, therefore, focuses on exploring the properties of these g-tori solutions, 
both analytically and numerically. In particular, we establish for g-tori energy localization 
laws analogous to those for g-breathers, using a semi-analytical approach. Our numerical 
experiments then show that such laws describe accurately the properties not only of exact 
g-tori solutions, but also of FPU trajectories with nearby initial conditions. 

Our work shares a common starting point with a recent paper by Giorgilli and Muraro 
where the authors have also explored the idea of FPU trajectories being confined on 
lower- dimensional manifolds embedded in the 2A^-dimensional FPU phase space. Their main 
result, proving that the confinement persists for times exponentially long in the inverse of 
e, is obtained in the spirit of the theory of Nekhoroshev [l^, [lol, using a variant of the 
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In 



formulation of the Nekhoroshev theorem for 'isochronous' systems 
fact, the theory of Nekhoroshev appears to offer quite a natural framework for studying 
analytically the metastability scenario. 



However, as has been pointed out quite early [25| a naive application of the Nekhoroshev 
theory in the FPU problem would break down as A^ ^ oo, since, under the assumption that 
the Nekhoroshev time T depends on the specific energy e, in the estimates T ~ exp(l/£:'^) of 
the theory the exponent is of the form c = 0{1/N) (see 2J] for a heuristic explanation). This 
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bad dependence of c on actually implies that T = 0(1) as ^ oo, hence Nekhoroshev's 
theory fails to predict long times for equipartition in that limit (a relevant theoretical result 
has only been obtained in lattices where a clear separation of the frequencies occurs into 



a low and high band; the low frequencies then become small parameters, see |26|). In 17 1 
one has instead c = 0(l/m), where m is one-half the dimension of the lower- dimensional 
manifold where confined FPU trajectories are expected to lie. This constitutes a signifi- 
cant improvement with respect to previous estimates, but still does not explain the natural 
packets correctly, since the latter's width was found in numerical experiments to vary as 



m oc e^^'^N 



131 ] ■ i.e. proportionally to A^, for fixed e. 



Our analytical theory relies on the use of the Poincare - Lindstedt method, through 
which we find scaling laws for the energy profile E{q) of a trajectory lying exactly on a q- 
torus. The consistency of the Poincare-Lindstedt construction on a Cantor set of perturbed 
frequencies (or amplitudes) is explicitly demonstrated. Numerically, we find that energy 
localization persists for appreciably long times, for trajectories neighboring a g-torus, even 
beyond the energy threshold where the g-torus becomes linearly unstable. The determination 
of linear stability for an s-dimensional torus (s > 1) is, of course, a subtle question, since no 
straightforward application of Floquet theory is available, as in the s = 1 case. Nevertheless, 
by employing an effective and reliable criterion for the stability of g-tori via the use of the 
recently developed method of 'Generalized Alignment Indices' (GALI, see 2!^), we are able 
to determine approximate critical parameter values at which a low- dimensional torus turns 
unstable, in the sense that orbits in its vicinity display chaotic behavior. 

It is important to remark at this point that, although the Poincare-Lindstedt approach 
is quite distinct from the Birkhoff method used in Nekhoroshev theory, it appears that the 
properties of exponential localization can be exploited to demonstrate (by subsequent use of 
the Birkhoff method), an even better behavior of the exponent c of Nekhoroshev estimates 
than can be found in the literature so far. A detailed exploration of this issue is deferred 
to a separate study. Nevertheless, a heuristic argument offered in the closing section of this 
paper suggests that the removal of the dependence of c on the number of degrees of freedom 
is possible, at least for trajectories close to g-breathers. 

Our paper is structured as follows: Section 2 presents our analytical results on the exis- 
tence and scaling laws of the 'g-tori'. We deal here only with the /5— case, but our approach 
can be readily extended to the a— case as well. We focus on the specific and most rele- 
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vant subset of 'g-tori' corresponding to zeroth order excitations of a set of adjacent modes 
(lo,i = i, i = I, ■ ■ ■ , s, whose energy profile E{q) is calculated analytically. Next, our an- 
alytical predictions are tested against numerical integration of specific orbits. Section 3 
examines the question of stability of g-tori and the persistence of energy localization of 
the FPU-trajectories when the linear stability of the 'underlying' g-tori is lost. Finally, we 
deal with the question of the long-term stability of exponentially localized FPU-trajectories, 
via heuristic estimates inspired by Nekhoroshev's theory. Section 4 summarizes the main 
conclusions of the present study. 



II. EXISTENCE AND STABILITY OF g-TORI 
A. The FPU /? model 

The /?— FPU Hamiltonian for a lattice of N particles reads: 

1 ^ 1 ^ /3 ^ 

k=l k=0 k=0 

where Xk is the k-th particle's displacement with respect to the equilibrium position and i/k 
is its canonically conjugate momentum. Fixed boundary conditions are defined by setting 

Xo = Xn+1 = 0. 

The normal mode canonical variables {Qq, Pq) are introduced by the linear canonical 
transformations 

q=l ^ ^ q=l ^ ' 



Substitution of ([2]) into ([T]) yields the Hamiltonian in the form H = H2 + Hi^ in which the 
quadratic part is diagonal 



H. = Y. ^^p^' (3) 

9=1 



with normal mode frequencies 



The quartic part of the Hamiltonian reads 

N 



(]g = 2sin(— 1^), l<g<iV . (4) 



Hi^ c^/j^ _|_ ^ ^ C q^i^jYi^ji^q^l^rn^nQ qQ iQ niQ n (5) 

q,l,m,n=l 
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where the coefficients Cq^i^rn,n take non-zero values only for particular combinations of the 
indices q, I, m, n, namely 

{lifg±/±m±n = 
-1 iiq±l±m±n = ±2{N + 1) 

in which all possible combinations of the ± signs are taken into account. Thus, in the new 
canonical variables, the equations of motion are: 

P ^ 

Qq ~l~ ^qQq 2fA^ ^ ^q,l,m,rS^q^l^in^nQlQmQn • (7) 

^ ^ l,m,n=l 

If /? = 0, the individual harmonic energies Eq = {P^ + Q'^Q'^)/2 are preserved by Eqs.([7]), 
i.e. the energies Eg form a set of integrals in involution. When /3 7^ 0, however, the 
harmonic energies become functions of time and only the total energy E = X]^i-^<?(^) 
conserved. The specific energy is then defined as e = E/N, while the average harmonic 
energy of each mode over a time interval < t < T is given by the integral Eq(T) = 
^j^Eg{t)dt. 

In classical FPU experiments, one starts with the total energy shared only by a small 
subset of modes. Then, for short time intervals T, we have Eq{T) ~ for all q corresponding 
to non-excited modes. Equipartition means that, due to the nonlinear terms, the energy 
will eventually be shared equally by all modes, i.e. 

lim Eq{T)=e , g = l,...,iV. (8) 

T — >oo 

The usual ergodic assumption of statistical mechanics leads to the conclusion that dHj) is 
violated only for orbits resulting from a zero measure set of initial conditions. The FPU 
paradox owes its name to the crucial observation that large deviations from the approximate 
equality Eg{T) ^ e occur for many other orbits as well. Depending on the initial conditions, 
these deviations are termed 'FPU recurrences' and are seen to persist even when T becomes 
very large. 



B. (7-tori and their construction by Poincare-Lindstedt series 



From the above discussion, one infers that the dynamics of FPU recurrences is governed 
by particular solutions of the FPU equations ([7]) for which energy localization occurs only 
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on a small subset of Fourier modes. Such solutions lie on tori of low-dimensionality, which 
we shall henceforth call 'g-tori', since they also turn out to be exponentially localized in 
Fourier space, like the 'g-breather' solutions discussed in l2|. 

We now introduce the main ingredients of our method of g-torus construction, using 
an explicitly solved example for = 8, whose solutions lie on a two-dimensional torus 
representing the continuation, for /3 7^ 0, of the quasiperiodic solution of the uncoupled 
{P = 0) system Qi{t) = Ai cosfiit, Q2 = A2CosQ2t, for a suitable choice of Ai and A2. 

To this end, we follow the Poincare - Lindstedt method and look for solutions Qq{t), 
q = 1, . . . ,8 expanded as series in the parameter a = P/2{N + 1), namely 

Q,it) = Q^^\t)+aQ^^\t)+a'Qf{t) + ..., g = l,...8 . (9) 

For the motion to be quasiperiodic on a two-torus, the functions Q^q\t) must, at any 
order r, be trigonometric polynomials involving only two frequencies (and their multiples). 
Furthermore, for the motion to represent a continuation of the unperturbed solutions Qi 
and Q2, the frequencies Ui and uj2 must be small corrections of the normal mode frequencies 
Qi, ^2- According to the Poincare - Lindstedt method, these new frequencies are also given 
by series in powers of a, as: 

ujq = Vlq + auj^^^ + a^J^^ + . . . g = l,2 . (10) 

The corrections are determined by the requirement that all terms in the differential equations 
of motion, giving rise to secular terms (of the form tsmuqt etc.) in the solutions Qq(t), be 
eliminated. 

Let us consider the equation of motion for the first mode, whose first few terms on the 
right hand side (r.h.s) are: 

Qi + nlQi = -aisntQl + enfnlQ^Ql + snln^QjQs + ...) (ii) 

Proceeding with the Poincare - Lindstedt series, the frequency Qi is substituted on the l.h.s. 
of Eq.lfTT]) by its equivalent expression obtained by squaring (JTOl) and solving for flf. Up to 
first order in cr this gives 

nj = - la^^J^^ + ... (12) 

Substituting the expansions ({9]) into (fTTj) . as well as the frequency expansion (fT2l) into the 
l.h.s. of (fTTj) . and grouping together terms of like orders, we find at zeroth order Q^^^ + 
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ujfQi = 0, while at first order 

- 3fi?fi3(Qf^)'Qf + ... (13) 

Repeating the above procedure for modes 2 and 3, we find that their zeroth order equations 
also take the harmonic oscillator form: 

qP+u;IqP = 0, QP+nlQP = . (14) 

Note that the corrected frequencies Ui, uj2 appear in the zeroth order equations for the 
modes 1 and 2, while the uncorrected frequency Vt^ appears in the zeroth order equation of 
the mode 3 (similarly, ^4, . . . , fig appear in the zeroth order equations of the modes 4 to 8). 
Continuing the construction of a solution which lies on a two-torus, we start from particular 
solutions of (with zero velocities at t = 0) which read: 

) (t) = Ai cos uit, Qf^ (t) = A2 cos uj2t, (t) = As cos figt 

where the amplitudes Ai,A2,As are arbitrary. If the solution is to lie on a two-torus with 
frequencies ui, U2, we must set As = 0, so that no third frequency is introduced in the 
solutions. In the same way, the zeroth order equations Q^g^ + ^gQ^q^ = for the remaining 
modes g = 4, ...8 yield solutions Q^q\t) = AgCosQqt and we set A4 = A5 = ... =^8 = 0. 
Thus, we are left with only two non-zero free amplitudes Ai, A2. 

Now, consider Eq. (fT3|) for the first order term Q^^\t). Only zeroth order terms Q^q\t) 
appear on its r.h.s., allowing for the solution to be found recursively. The crucial remark is 
that by the choice A3 = . . . = Ag = 0, one also has Qf\t) = . . . = Qf\t) = 0, whence only 
a small subset of the terms appearing in the original equations of motion survive on the r.h.s. 
of Eq. (|T3l) . namely those in which none of the functions Qf\t), . . . , Qf \t) appears. As a 
result, equation (fT3|) is simplified dramatically and upon substitution of Qf\t) = Ai coscuit, 
Q'^\t) = A2Cosuj2t reduces to: 

^ + LulQ^^'^ = 2fiicuf ^ cos tuit - Sn^Al cos^ uit 

— GQlQlAiAlcosuit cos^ Lj2t ■ (15) 

This can now be used to fix 00^^ so that no secular terms appear in the solution, yielding 

= IaIqI + ^Aln.Ql 
9 



while, after some simple operations, we find for Q^^^: 



(1) _ 3A?0f cos Scjit 3AiAlnlQlcos{uJi + 2uj2)t 
~ 3%Jl ^ 2[{uJi + 2uj2y -uf] 

SAiAlnjnl cos(u;i - 2uj2)t 
+ 2[{uj,-2u2y-ujf] ■ ^ 

By the same analysis, we fix the frequency correction of the second mode: 

4'^ = -AjQl + -Alnln2 
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and obtain the solution 



(1) _ 3A3fi4 cos Scjat SAfAsfifn^ cos(2u;i + a;2)t 
SA^Aa^i^^ cos(2u;i - cja)^ 



(17) 



2[{2uj^-uJ2y-uji] 

which has a similar structure as the first order solution of the first mode. For the third order 
term there is no frequency correction, and we find 

(1)^^^ A^Q^Q^ fScosuJit ^ cos3c(Jit 



4 9cj2-^]2 



3741742^1^2^3/ cos(u;i — 2a;2)t cos (c^i + 2^012 )i 2cosa;it\ , . 

^ 4 V K - 2cu2)2 - ^1 ^ K + 2cu2)2 -Qj^ Ljf -Ql) ■ ^ ' 

We may thus proceed to the sixth mode to find solutions in which all the functions 
Qf\ . . . , Qf^ are equal to zero, while the functions Q^^\ . . . , Q^^^ are non zero. However, a 
new situation appears when we arrive at the seventh and eighth modes. A careful inspection 
of the equation for the term Q'^\t) 

8 

Q«+fi2Q«=- J2 C-j^l^^^nWl^AQf^Qt^Qf (19) 
i,m,n=l 

shows that there can be no term on the r.h.s. which does not involve some of the functions 
Qf \ . . . , Qf^- Again, this follows from the selection rules for the coefficients of Eq.([H]). Since 
all these functions are equal to zero, the r.h.s. of Eq. flT^ is equal to zero. Taking this into 
account, we set Q^^\t) = 0, so as not to introduce a third frequency ^7 in the solutions, 
whence the series expansion (|9]) for Qi{t) necessarily starts with terms of order at least 
O(cr^). The same holds true for the equation determining Q^^\t). 
Some remarks regarding the above construction are in order: 
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i) Consistency. The solutions (fT6|) - (fTSll (and those of subsequent orders) are meaningful 
only when the frequencies appearing in the denominators satisfy no commensurability con- 
dition. The spectrum of uncorrected frequencies fig, given by (jl]), is fully incommensurable 
only if is either a prime number minus one, or a power of two minus one (285. In all 
other cases, there are commensurabilities among the unperturbed frequencies, examples of 
which are given e.g. in 29|]. However, such commensurabilities do not affect the consistency 
of the construction of the Poincare-Lindstedt series, because it can be shown that no di- 
visors of the form "^q^q appear in the series at any order and for any integer vector 
n = (ni, n2, . . . , Uq) 7^ (0, 0, . . . , 0). The proof of this statement follows from the fact that 
the kernel differential equations (like Eq. (IT^ ) determining all terms Q^q \ q = 1,...8 at 
order k read: 

fc-1 

cos^,t+ Yl 41,n2 cos[(ni^i + n2a;2)t] ifg = l,2 (20) 

l"ll + l"2l#0 

while for g = 3, ...,8 the same equation holds with B^'' = and uj"^ on the l.h.s. replaced 
by fig. The coefficients Bq''^ and Ag^2i,n2 in (120|) are determined at previous steps of the 
construction (see Eq. fl2^ below for arbitrary s). It follows, therefore, that all new divisors 
appearing at successive orders, belong to one of the following sets: 

(ni ± l)uji + n2Co'2, niui + (n2 ± 1)^2, fig ± {niUJi + n2CJ2), g = 3, . . . 8 

for ^1,^2 G Z, l^il + \n2\ 7^ 0, whence we deduce that no zero divisors can ever show 
up in the Poincare-Lindstedt series due to commensurabilities between the unperturbed 
frequencies fig. On the other hand, we can also exclude the appearance of zero divisors due 
to the perturbed frequencies uji,uj2, if the frequencies in a quasiperiodic solution are fixed in 
advance so that the two frequencies do not belong to the countable set of all planes defined 
by the relations (rii ± l)^^;! + r;,2Ci;2 = 0, niUi + (n2 ± l)ci;2 = 0, fig + {niUi + n2UJ2) = 0, 
q = 3,... 8, for all integer values ni,n2 G Z, \ni\ + \n2\ 7^ 0. We stress that fixing the 
frequencies in advance is a necessary ingredient of the Poincare - Lindstedt method, since 
otherwise the equations at all orders (as, for example, Eqs. (|T6l) . (|T71) and (|T8l) ) would not be 
solvable. Furthermore, this frequency specification is analogous to the procedure followed in 
the construction of Kolmogorov's normal form representing solutions on KAM tori (see 
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for a detailed comparison of the two methods). Since the frequencies ui, uj^ are functions 
of the amphtudes Ai,A2, one deduces that the formal consistency of the method can be 
established in the complement of all excluded planes, i.e. a Cantor set in either the frequency 
space (ti;i,ci;2) or the amplitude space {Ai,A2). 

The above demonstration of consistency is readily generalized to the construction of s- 
dimensional g-tori by the nonlinear continuation of the set of linear modes qi, i = 1, ...,s. 
Namely, one can demonstrate that the consistency holds on a Cantor set of amplitude 
multiplets (Ag^, Ag^, . . . , AgJ such that the resulting perturbed frequencies ujg-,i = 1, . . . , s 
not lie on any one of the planes: 

(ni + mi)ujci^ + (n2 + m2)ujg2 + . . . + (n^ + m^)t^g^ = 

or 

rig ± {niUJg^ + 722^92 + • • • + nsUJqs) = 

where g = {1, . . . , A^}\{gi, . . . ,qs}, ^i, . . . , G A^, |r;,i| + . . . + \ns\ ^ and rUg = -1, 
or 1 for all g = 1, . . . , s. For example, if the first s modes are excited by amplitudes A^, 
k = 1, ...s, the formula for the perturbed frequencies reads 

Ug = ng + ^ngJ2Al^l-lAK + 0{cr'-A,,...A,), q = l,...,s (21) 

k=l 

Fixing the values of the frequencies Ug in advance implies that Eqs.( l2T|l should be regarded 
as yielding the (unknown) amplitudes A^ for which the quasiperiodic solution exhibits the 
chosen set of frequencies on a Cantor set in the space of the amplitudes Ai, . . . ,As. The case 
s = 1 implies that the same property holds for the Poincare - Lindstedt series representing 
g-breathers 21]. That is, while by Lyapunov's theorem the continuation of the periodic 
orbits is guaranteed in an open domain of values AgQ, Poincare-Lindstedt series can only be 
constructed on a Cantor subset of this domain. Yet, this is sufficient for our present purpose, 
which is to determine by a semi-analytical approach scaling laws for the energy localization 
on g-tori (or g-breathers). 

ii) Convergence. Even after consistency is demonstrated, no proof has yet been pro- 
vided for the convergence of the series. As demonstrated in the works of Eliasson [31] and 



Gallavotti [32, |33|], the question of convergence of the Lindstedt series is notoriously diffi- 



cult even in simple Hamiltonian systems. This is because the Lindstedt series for orbits on 
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FIG. 1: Comparison of numerical (points) versus analytical (solid line) solutions, using the Poincare 
- Lindstedt series up to order O(fT^), for the temporal evolution of the modes (a) q = 1, (b) q = 3, 
(c) q = 7, when Ai = 1, A2 = 0.5, and N = 8, (3 = 0.1 and (d) the time evolution of the indices 
GALI2 to GALIq up to a time t = 10^ show that the motion lies an 2-dimensional torus. 



invariant tori are convergent, but not absolutely (see the review by Giorgilli 3J]). On the 
other hand, it is possible to make an absolutely convergent classical expansion by use of the 
Kolmogorov normal form as developed by Giorgilli and Locatelli 35|]. Such analysis, how- 
ever, is quite cumbersome and will be deferred to another publication as it would obscure 
the results presented here. 

Thus, we prefer to justify our previous statements by numerical simulations, taking our 
initial conditions from the analytical solutions ([9]) at t = and using the GALl method 27 1 
to demonstrate that the solutions lie, indeed, on two-dimensional tori, as shown in Figure [H 
The numerical solution for the modes Qi{t), Qsit) and Qjit) is checked against the analytical 
solution via the Poincare-Lindstedt series, truncated at second order with respect to a, when 
N = 8, jS = 0.1, and Ai = 1, A2 = 0.5. The size of the error is found to be precisely that 
expected by the truncation order in FigH^-c, while in FigUJi the GALI method shows that 
the numerical orbit lies on a 2-torus. 
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Let us recall that according to Skokos et al. [271], the indicators GALI^, k = 2,3,..., 
for chaotic orbits, decay exponentially fast, due to the attraction of all deviation vectors 
by the most unstable direction corresponding to the maximal Lyapunov exponent. On the 
other hand, if an orbit lies on a stable s-dimensional torus, the indices GALI2, . . . , GALIg 
oscillate about a non-zero value, while the GALIg^j, j = 1, 2, .. follow asymptotically power 
laws falling at least as . 

This is precisely what we observe in Fig. [T]i. Namely, after a transient initial interval 
(required for phase mixing to become effective), the GALI2 index stabilizes at a constant 
value GALI2 — 0.1, while all subsequent indices, starting from GALI3 decay following a 
power law as predicted by the theory. The time for phase mixing is estimated to be of order 
N^/j3 ~ 10^. This implies that beyond t = 10^ we should observe the expected asymptotic 
behavior of the GALI indices to set in, as is the case in Fig{T]i. Thus, we conclude that 
the motion lies on a 2-torus, exactly as predicted by the Poincare-Lindstedt construction, 
despite the fact that some excitation was provided initially to all modes. 

(iii) Presence and accumulation of small divisors. There are, of course, small divisors 
appearing in terms of all orders beyond the zeroth. First, the low-mode frequencies satisfy 
ujq ~ TTq/N, and hence divisors like ujf, or {ui — 2002)'^ ± ujf (appearing e.g. in Eq. f|T6l) ) are 
small and care must be taken regarding their effect on the growth of terms of the series 
at successive orders. In fact, the most important effects are introduced by nearly resonant 
divisors, like dujf — in Eq. lfTSl) . Since the first order corrections of the frequencies ooi and 
UJ2 are 0{(3A'j/N'^), for Aj, (3 sufficiently small, one may still use for them the approximation 
given by the first two terms in the sinus expansion of (jll) namely: 

- (iV + 1) ~ 24(iV + l)3 ^^^^ 

the error being 0{A^(3/N^) for uji,uj2, and 0{{q/N)^) for all other frequencies. This implies 
that a divisor like 9uf — Q"^ can be approximated by the relation 

.2, ,2 02, oM., 2g7r TT^iq^ -q) ttY , f ^l' 



1 1^"^ 1^ 1 gj\ (iv + 1) 24(Ar + 1)3 12m \m 

In general, the terms produced at consecutive orders involve products of divisors, whose 
influence on the size of the series terms must be taken into account as regards estimates 
of the profile of energy localization for the g-tori solutions, as explained in subsection II. D 
below. 
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FIG. 2: Same as in Fig. [H showing the existence of a 4-torus in the system with = 16, (3 = 0.1, 
when Ai = 1, A2 = 0.5, ^3 = 0.333..., A4 = 0.25. The temporal evolution Qq{t) is shown for the 
modes {a)q = 1, (b)g = 5, {c)q = 13. (d) The GALI^, /c = 2,3,4 indices, are seen to stabilize after 
t ~ 10^, while for k > 5 the indices continue to decrease by power laws. 



(iv) Sequence of mode excitations. The profile of energy localization along a g-torus 
solution is determined by the sequence of mode excitations arising as the recursive scheme 
proceeds to subsequent orders. The term 'excitation' here means that the solutions of the 
Poincare - Lindstedt method should be non-zero for the first time at the order where it is 
claimed that the excitation takes place. For example, as already explained in our previous 
construction of a 2-torus solution, the modes 1,2 are excited at zeroth order, the modes 3, 4, 
5, and 6 at first order, and the modes 7 and 8 at the second order of the recursive scheme. 

Figure [2] presents a one more example showing the comparison between our analytical and 
numerical results for the modes Qi{t), Q^it), and Quit), along a 4-torus solution constructed 
precisely as described above, with N = 16, (3 = 0.1 and by exciting modes 1 - 4 at zeroth 
order, via the amplitudes Ai = 1, A2 = 0.5, = 0.333 . . ., ^4 = 0.25. In this case, we find 
that the modes excited at the first order of the recursive scheme are g = 5 to g = 12, while 
the modes excited at second order are g = 13 to g = 16 and the GALI method shows that 
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the motion occurs on a 4-torus, since GALI^ is the first index to drop asymptotically like 
t-^ (see Fig. EH). 

The sequence of excitation of different modes plays demonstrably a crucial role in esti- 
mating the profile of energy localization, since the amplitudes of all excited modes at the 
r-th order have a pre-factor = {(3/2{N + 1))'^. This is the subject of the next subsection, 
in which a proposition is provided for the sequence of mode excitations, in the generic case 
of arbitrary N and arbitrary dimension s of the low-dimensional torus (with the restriction 
s < N). Next, the consequences of this proposition are examined on the localization profile 
of g-tori and nearby FPU trajectories. 

C. Sequence of mode excitations 

In order to motivate the results of this subsection, let us first limit ourselves to what 
happens in the case of the solutions of Figures [Hand El Figure |3] shows the average harmonic 
energy of each mode over a time span T = 10^ in the cases of the g-torus of Fig. [1] and Fig. 
El shown in Fig. [3^ and Fig. [3)d respectively. The numerical result (open circles) compares 
excellently with the analytical result (stars) obtained via the Poincare - Lindstedt method. 
The filled circles in each plot represent 'piecewise' estimates of the localization profile in 
groups of modes excited at consecutive orders of the recursive scheme. The derivation 
and exact meaning of these theoretical estimates will be analyzed in detail below. Here 
we point out their main feature, showing a clear-cut separation of the modes into groups 
following essentially the sequence of excitations predicted by the Poincare-Lindstedt series 
construction. Namely, in Fig. [3^ we see clearly that the decrease of the average energy 
Eg along the profile occurs by abrupt steps, with three groups formed by nearby energies, 
namely the group of modes 1,2, then 3 to 6, and then 7,8. The same phenomenon is also 
seen in Fig. [3b, where the grouping of the different modes follows precisely the sequence of 
excitations (1 to 4), (5 to 12), (13 to 16) as predicted by the theory. 

Let us, therefore, extend our approach to a more general consideration of the structure 
of the solutions lying on low-dimensional tori. Using the Poincare - Lindstedt method, a 
solution and its frequencies are expanded in series of the form 

oo oo 

Q,{t) =J2^'Qi'\t), = Y^a'uf^ (23) 

fc=0 fc=0 
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FIG. 3: (Color online) The average harmonic energy Eg of the q-th mode as a function of q, after a 
time T = 10^ for (a) the 2-torus and (b) 4-torus solutions (open circles) corresponding to the initial 
conditions used in Figs|T] and [2] respectively. The stars are Eg values calculated via the analytical 
representation of the solutions Qg{t) by the Poincare -Lindstedt series. The filled circles show a 
theoretical estimate based on the average energy of suitably defined groups of modes (see Eg. ([35]) 
and relevant discussion in the text). 

where Uq^^ = fig. Substituting fl23l) into the equations of motion ([7]) and separating terms 
of like orders, the equations at order k read 

k n\ 

Qf+ulQf^ = 5^(^ a;(-)a;("))gf--) - 

"1=1 n2=0 
N k-1 

Qg mm^nCqlmn J] Q^"^ Qt'^ Q^n'^ ■ (24) 

/,m,n=l "-1,2,3=0 

n-i +71-2 +T13— fc — 1 

Note that Eqs.f l24p are still quite general. Let us consider, therefore, the case where 
only a subset of modes 1 < qi < q2 < ■■■ < qs < N are excited at the zeroth order of the 
perturbation theory, assuming that Q^q^ ^ iff g G {qiiq2-, ■ ■ ■ iQs]- The modes gi, . . . , 
need not be consecutive. We wish to see how this type of zeroth order excitation propagates 
at subsequent orders. More specifically, we wish to determine for which q values one has 
Q{k <k) _ Q{k) _^ modes q are first excited at the k-th order. The answer is 

provided by the following 



Proposition: Let the starting terms of a Poincare-Lindstedt series solution with s 
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frequencies be set as 

Qf} = Ag^cosiugi + 4>g^) /o F i = 1 , 2 , . . . , s , 1 < ^1 < ^2 < • • • < < 

Qf = for all q ^ . (25) 

Then, besides the terms Q^_{t), the Poincare-Lindstedt series terms Q^q\t) which are per- 
mitted to be non-zero at the kth order of the series expansion are given by the values of 
q = q^^^ satisfying 

q^^'^ =\ 2X{N + 1) - ruk \, mk{mod{N + 1)) ^ (26) 

where ruk can take any of the values \ ±qi-^ ± qi^ ± ... ± qi^^^-^ \, with ii,i2, ...,i2k+i ^ {1, s} 
for any possible combination of the ± signs, and A = [{nik + N)/2{N + 1)]. 



The proof of the proposition is given in Appendix A. Some simple examples clarify 
the use of the rule (|26l) : 

q -Breathers: If we excite only one mode qi at zeroth order, new modes are excited one 
by one at subsequent orders and one has = {2k + 1)^1. As long as < A^, the newly 
excited modes are g*^'^^ = mk = {2k + l)gi, i.e., exactly as predicted by Flach et al. A 
particular case arises when qi = 2{N + l)/3, (A^ + l)/2 or (A^+ l)/3. Then, one readily sees 
that no new modes are excited at any subsequent order, which is in agreement with a well 



known result 



36l | (see also j37|). 



2-dimensional q-tori: Assume we excite the modes gi = 1, ^2 = 2 at zeroth order. At 
first order {k = 1) we have again q = ± qi^ ± qi^\ with 21,22, "^3 € {1;2}. We readily find 
that the newly excited modes are g = l + l + l = 3, g = l + l + 2 = 4, g=l + 2 + 2 = 5, and 
g = 2 + 2 + 2 = 6. At order k = 2, the first newly excited mode isg = l + l + l + 2 + 2 = 7, 
while the last newly excited mode isg = 2 + 2 + 2 + 2 + 2 = 10. In general, at order k > 1 
the newly excited modes are 2{2k — 1) + 1 < q < 2{2k + 1). 

s -dimensional q-tori: Assume we excite the modes gi = 1, ^2 = 2, . . ., = s at the 
zeroth order. In the same way as above we find that the newly excited modes at order k > 1 
are s{2k - 1) + 1 < g < s{2k + 1). 

More complicated choices of the initially excited modes gi, . . . , g^ lead to very interesting 
localization patterns that will be the subject of a separate study. 
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D. Profile of the energy localization 



In order to study now in Fourier space energy localization phenomena associated with 
FPU-trajectories, let us construct estimates for all Qq{t) terms participating in a particular 
g-torus solution in which the s first modes (gi, ^2, • • • , Is) = (1; 2, . . . , s) are excited at zeroth 
order of the theory, with amplitudes Ai, A2, . . . , As respectively. 

Denoting by q^''^ the indices of all modes which are newly excited at A;-th order, according 
to the proposition of the previous subsection, one has 

e = {max(l, (2A; - l)s + 1), . . . , (2A; + , (27) 

whence the following useful estimates are obtained 

Vg(^) G Ms,,, k > 1 g('=) < {2k + l)s, < ^^^^^ ■ (28) 

Let us now define the 'majorant' norm in the space of trigonometric polynomials / 

11/11 = El^'^l (29) 

k 

where A^ are the coefficients of a trigonometric polynomial with, say (for simplicity), only 
cosine terms 

/ = 5^AkCos(k-0) . (30) 

k 

The quantity ( !29l) satisfies all properties of the norm. Since the series terms of a mode 
Qq(k){t) satisfy Q^^l) = 0, Vra < k, Q^^d) 7^ 0, from Eq.(|2^ we deduce that the equation 
determining reads 

n(^) _L , ,2 n^'^) _ 

.-1 / \ 



(31) 



EO O O /n("i) /n("2) n^n-i) 
i ',("2) i ^q("3) V^(ni) V^(n2) Vq{n3) 

"1,2,3=0 \ (g("l),g{".2),q("3))GD J 

where 

T^q = { iljl , (Ih , 9j3 ) : ? ± ± ?J2 ± % = 0} . 

Assuming that the Lindstedt series are convergent, omission of explicit reference to higher 
order terms Q^l^(t)^\ Q^l^(t^\ • • •, is justified as long as only an estimate of the size of the 
oscillation amplitude of the mode g'^'^^ is sought. 
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The average size of the oscillation amplitudes of each mode along an s-dimensional q- 
torus follows now from estimates on the norms of the various terms appearing in Eq. flST]) . If 
we denote by A^'^^ the mean value of all the norms ||(5„(fe)||, we find the following estimate: 



'9' 



(32) 

where Aq = A^^'^ is the mean amplitude of the oscillations of all modes excited at the zeroth 
order of the perturbation theory, and C is a constant of order 0(1). The proof of (132!) is 
deferred to Appendix B, in which the analytical estimate C ^ 3/2 isgiven. In fact, fl32l) is a 
straightforward generalization of the estimate given in Flach et al. [2] for g-breathers, while 
the two estimates become identical (except for the precise value of C) if one sets s = 1 in 
Eq. (l32i) . and go = 1 in Flach's g-breathers formulae. 

Thus, the g-tori provide an explanation for the results reported in [l^ and 11 1 and offer 



a bridge between the 'natural packet' approach and the interpretation of energy localization 
for FPU-trajectories based on g-breathers. In particular, the physical picture suggested by 
the above analysis is that, starting with initial conditions near g-breathers, a 'backbone' is 
formed in the phase space by a hierarchical set of solutions which are, precisely, the solutions 
lying on low- dimensional g-tori (of dimension s = 1,2, ...,s << A^). All FPU-trajectories 
with initial conditions within this set exhibit a profile of the energy localization characterized 
by a 'stepwise' exponential decay, with step size equal to 2s, as implied by Eq. (!27|) . All the 
modes g^'^-' E ■Mk,s share a nearly equal mean amplitude of oscillations, which follows the 
estimate of Eq. (|32il . 

Using ( l32l) . we find it convenient to obtain 'piecewise' estimates of the energy of each 
group using a formula for the average harmonic energies E^''^ of the modes g^'^-'. To achieve 
this, note that the total energy E given to the system can be estimated as the sum of the 
energies of the modes 1, . . . , s (the remaining modes yield only small corrections to the total 
energy), i.e. 



E ~ SUJimAn ~ 



On the other hand, the energy of each mode g^'^-' can be estimated from 



2 ''''' \2{N + 1) J ' ' 22'=+i(Ar + l)2fc+2 
which, in terms of the total energy E^ yields 



^„^^(CWJ^)' . (33) 
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Once again, the similarity of Eq. (l33|) with the corresponding equation for g-breathers is 
obvious. The latter equation reads [2I 

i?(2.+i),o - E,, ^ j (34) 

where go is the unique mode excited at zeroth order of the perturbation theory. Note, in 
particular, that the integer s plays in Eq. (!33l) a role similar to that of go in Eq.( !34l) . This 
means that the energy profile of a g-breather with go = s presents the same exponential law 
as the energy profile of the s-dimensional g-torus. But the most important feature of the 
latter type of solutions is that the profile remains unaltered as increases, provided that: 
i) a constant fraction M = s/N of the spectrum is initially excited, (i.e. that s increases 
proportionally to A^), and ii) the specific energy e = E/N remains constant. Indeed, in 
terms of the specific energy e, ( !33l) takes the form 

Ei^>^±(^r (35) 

i.e. the profile becomes independent of A^. A similar behavior is recovered in the g-breather 
solutions provided that the 'seed' mode go varies linearly with A^, as was shown in detail in 
refs.yisg]. 



E. Numerical results for FPU-trajectories 

Examples of the 'stepwise' profiles predicted by Eq. fl33|) in the case of exact g-tori solu- 
tions are shown by filled circles in Fig. [3l concerning the solutions depicted in Figs. [1] and [2] 
(in all fittings we set C = 1 for simplicity). From these one can see that the theoretical 'piece- 
wise' profiles yield nearly the same average exponential slope as the profiles obtained either 
numerically, or analytically by the construction of the solutions via the Poincare - Lindstedt 
method. Thus, the estimates (|33|) . or fl35l) . appear quite satisfactory for characterizing the 
localization profiles of exact g-tori solutions. 

The key question now, regarding the relevance of the g-tori solutions for the interpretation 
of the FPU paradox, is whether Eqs. (|33l) or (|35l) retain their predictive power in the case 
of generic FPU trajectories which, by definition, are trajectories started close to, but not 
exactly on a g-torus. 

An answer to this question is partly contained in the results presented in Figures H] and 
|5l Figure H] shows the energy localization profile in numerical experiments in which (3 is 
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FIG. 4: (Color online) The average harmonic energy Eg of the g-th mode over a time span T = 10^ 
as a function of q in various examples of FPU-trajectories, for /3 = 0.3, in which the s(= A^/16) 
first modes are only excited initially via Qq{0) = Ag, Qq{0) = 0, q = 1, . . . , s, with the Ag selected 
so that the total energy is equal to the value E = H indicated in each panel. We thus have (a) 
AT = 64, ^ = 10-^ (b) N = 128, E = 2x IQ-^ (c) N = 256, ^ = 4 x IQ-^, (d) N = E = IQ-^, 
(e) N = 128, E = 2x IQ-'-^, (f) N = 256, E = i x 10~^. The specific energy is constant in each of 
the two rows, i.e. e = 1.5625 x 10^^ in the top row and e = 1.5625 x 10^^ in the bottom row. The 
dashed lines represent the average exponential profile Eg obtained theoretically by the hypothesis 
that the depicted FPU trajectories lie close to g-tori governed by the profile (j35|) . 

kept fixed {(3 = 0.3), while takes the values N = 64, N = 128 and = 256 (although in 
subsection IIB the derivation of explicit g-tori solutions was practically feasible by computer 
algebra only up to a rather small value of (A^ = 16), in the present subsection the results 
with numerical trajectories are extended to much higher values of A^). In all six panels of 
Fig. mthe FPU-trajectories are computed starting with initial conditions in which only the 
s = 4 (for A^ = 64), s = 8 (for A^ = 128) and s = 16 (for A^ = 256) first modes are excited at 
t = 0, with the excitation amplitudes being compatible with the values of the total energy 
E indicated in each panel, and constant specific energy e = 1.5625 x 10~^ in the top row 
and e = 1.5625 x 10~^ in the bottom row of Fig. HI 
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The resulting trajectories differ from g-tori solutions as follows: In the g-tori all modes 
have an initial excitation, whose size was estimated in Eqs. fl32|) or fl33l) . while in the case of 
the FPU trajectories only the s first modes are excited initially, and one has Qq{0) = for all 
modes q > s. Furthermore, since in the g-tori solutions one also has \ \Qq>sit) 1 1 << \ \Qq<s(t) 1 1 
for all t, the FPU trajectories can be considered as lying in the neighborhood of the g-torus 
solutions, at least initially. The numerical evidence is that if E is small, they remain close 
to the g-tori even after relatively long times, e.g. t = 10®. 

This is exemplified in Fig. IH in which one sees that the average energy profiles of the 
FPU-trajectories (at t = 10®) exhibit the same behavior as predicted by Eg. (1331) . for an 
exact g-torus solution with the same total energy as the FPU trajectory in each panel. For 
example, based on the values of their average harmonic energy, the modes in Fig. (in 
which s = 4) are clearly separated in groups, (1 to 4), (5 to 12) and (13 to 20), etc., as 
foreseen by Eg. (1271) for an exact 4-torus solution. The energies of the modes in each group 
have a sigmoid variation around a level value characteristic of the group, which is nearly the 
value predicted by Eg.( |33l) . The grouping of the modes is distinguishable in all the panels of 
Fig. m a careful inspection of which verifies that the grouping follows the laws found above 
for g-tori. Also, if we superpose the numerical data of the three top (or bottom) panels we 
find that the average exponential slope is nearly identical in all the panels of each row, a 
fact consistent with Eg. fl35|) . according to which, for a given fraction M of initially excited 
modes, this slope depends on the specific energy only, i.e. is independent of for constant 
e. 

When we increase the energy, the FPU-trajectories resulting from s initially excited 
modes start deviating from their associated exact g-tori solutions. As a conseguence, the 
energy profiles of the FPU-trajectories start also deviating from the energy profiles of the 
exact s-tori. This is evidenced by the fact that the profiles of the FPU-trajectories become 
smoother, and the groups of modes less distinct, while retaining the average exponential 
slope predicted by Eg. (l35l) . This 'smoothing' of the profiles is discernible in Fig. [5K, in 
which the energy is increased by a factor 50 with respect to Fig. HJi, for the same values of 
A^ and (3. Also, in Fig. [5^ we observe the formation of the so-called 'tail', i.e. an overall rise 
of the localization profile at the high-freguency part of the spectrum, accompanied by spikes 
at particular modes. This is a precursor of the evolution of the system towards eguipartition, 
which manifests itself earlier in time for larger energies. 
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FIG. 5: (Color online) Same as in FigS^ but for larger energies, namely (a) E = 0.05, (b) E = 0.1, 
(c) E = 0.2, (d) E = 0.3, (e) E = 0.4, (f) E = 0.5. Beyond the threshold E ~ 0.05, theoretical 
profiles of the form (j33p yield the correct exponential slope if s is gradually increased from s = 4 
in (a) and (b) to s = 6 in (c), s = 7 in (d), s = 10 in (e), and s = 12 in (f). 

Nevertheless, the important remark is that the phenomenon of exponential localization 
of the FPU trajectories persists, and is still characterized by laws like Eq.( l33l) . even when 
the energy is substantially increased. Furthermore, at energies beyond a threshold value, an 
interesting phenomenon occurs which is worth mentioning: For fixed (see e.g. Figure [5l 
where N = 64), as the energy increases, a progressively higher value of s needs to be used 
in Eq. fl33|) . so that the theoretical profile yields an exponential slope that agrees with the 
numerical data. 

As evidenced in Fig. [5] for f3 = 0.3, = 64, the threshold is E ^ 0.1. This value splits 
the system in two distinct regimes: One for E < 0.1, where the numerical data are well 
fitted by a constant choice of s = 4 in Eq. (!33|) (indicating that the FPU-trajectories are 
indeed close to 4-tori), and another for E > 0.1, where 'best-fit' models of Eq. fl33|) occur for 
values of s increasing with the energy, e.g. s = 6 for E = 0.2, rising to s = 12 for E = 0.5. 
This indicates that the respective FPU-trajectories are close to g-tori with a progressively 
higher value of s (with s > 4), despite the fact that only the four first modes are excited by 
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the initial conditions of these trajectories. 

This behavior is analogous to the 'natural packet' scenario described by Berchialla et al. 
[13I . in which a set of modes is seen to share the energy after some time even if this energy 
is initially given to only the first mode. These authors also observed that the law giving 
the localization profile of their metastable states stabilizes as the energy increases. Indeed, 
according to Eg. (1331) . such a stabilization implies that in the second regime the width s 
depends asymptotically on E as s (x E^^"^, or, from Eq. fl35|) that M oc e^^^. This agrees well 
with estimates on the width of natural packets formed by the /3— FPU model described in 
12|. 



III. STABILITY OF THE MOTION NEAR g-TORI 



A. Linear stability 

The linear stability of g-breathers can be studied by the implementation of Floquet theory 



see 



m 



which demonstrates that a g-breather is linearly stable as long as 

R^^JME±21^r^0(l/N^) (36) 
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This result is obtained by analyzing the eigenvalues of the monodromy matrix of the lin- 
earized equations about a g-breather solution constructed by the Poincare - Lindstedt series. 
Numerical verification can also be used to analyze the dynamics about the fixed point cor- 
responding to a g-breather under the Poincare map of the fiow of Eqs.Q. 

In the case of g-tori the above techniques are no longer available. Nevertheless, a re- 
liable numerical criterion for the stability of g-tori is provided by the GALI indices j27|. 
According to this method, if a g-torus becomes unstable beyond a critical energy threshold, 
the deviation vectors of trajectories started exactly on the g-torus are attracted by its un- 
stable manifold, whence all the GALI indices beyond and including GALI^^i {k being the 
dimension of the unstable manifold) fall exponentially fast. Naturally, if one starts with tra- 
jectories in the vicinity of an unstable g-torus, these trajectories are weakly chaotic. Thus, 
it turns out that all GALI indices start falling exponentially after a transient time, and this 
can be checked by calculating the time evolution of the lowest index, i.e. GALI2. 

Using this criterion, we examined the stability of g-tori and determined approximately 
the value of the critical energy E^ at which a g-torus turns from stable to unstable. Figure 
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FIG. 6: (Color online) Assuming that the critical energy Ec, at which the GALI2 index shows that 
a (7-torus destabilizes is fitted by the law Ec = Al3~^, we show in the upper curve (triangles) the 
dependence of ^ on for an FPU trajectory started by exciting initially only the q = 1,2 modes. 
The middle curve (filled circles) corresponds to a similar calculation for FPU trajectories started 
near a q'-breather solution, where only the q = 1 mode is excited. The dashed line corresponds to 
A ~ N^^, according to the law of Eg. ([55]) . 

|6] shows an example of this calculation using FPU trajectories started close to a 2-torus. 
The quantity A shown in the ordinate corresponds to a calculation keeping N fixed and 
varying /3, until a critical energy Ec is determined, beyond which the GALI2 index loses 
its asymptotically constant behavior. All calculations refer to a maximum time t^ax = 10^ 
up to which we require that the exponential fall-of of GALI2 must have been observed (in 
general Ec decreases as tmax increases, tending to an asymptotic limit as tmax 00). The 
values of Ec found this way provide upper estimates for the transition energy at which the 
exact g-torus turns from stable to unstable, i.e., for any choice of tmax the transition energy 
is lower than the value of Ec found by the GALI2 method. 

As expected, due to the obvious scaling of the FPU Hamiltonian by (3, the critical energy 
for all considered values of turns out to be well fitted by a power law Ec = Af3~^. However, 
the fitting constant A depends also on A^, as seen in Fig. [61 From the triangle data in the 
figure it is clear that the dependence of A on A^ is weaker than A^~^, i.e. the law predicted 
by (|36|) for the g-breathers. 
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FIG. 7: (Color online) (a) Time evolution of the GALI2 index up to t = 10'' for an FPU trajectory 
started by exciting the q = 1 and q = 2 modes in the system N = 32, (3 = 0.1, with total energy 
E = 2, i.e. higher than Ec = 1.6. (b) Instantaneous localization profile of the FPU trajectory of 
(a) at t = 10^ . (c) Same as in (a) but for N = 64, E = 1.323 (in this case the critical energy is 
Ec = 1.24). (d) Same as in (b) but for the trajectory of (c). 

On the other hand, the dependence A oc A^~^ was approximately found for FPU- 
trajectories started close to a g-breather, where only the mode g = 1 is initially excited. In 
that case the numerical data (filled circles) have a slope close to that predicted by Eq. (!36!l . 
although the whole numerical curve is shifted upwards with respect to the dashed line, a 
fact verifying that the GALI2 method yields critical energies E^. which are higher than the 
value at which the g-breather becomes unstable. 

The weak dependence of A on on the upper curve of Fig. [6] is a numerical indication 
that the g-tori solutions are more robust than the g-breathers regarding their linear stability. 
This numerical behavior is related to the results of the previous subsection, which indicate 
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that the number of modes excited initially does not always coincide with the dimension of 
the g-torus which the FPU-trajectory eventually approaches. Also, the destabilization of 
simple periodic orbit represented by a g-breather does not imply that the tori surrounding 
the breather are also be unstable. 

At any rate, the most important remark concerning the linear stability of g-breathers or q- 
tori is that the exponential localization of FPU trajectories persists even after the associated 
g-breathers or g-tori have been identified as linearly unstable by the GALI criterion. This 
behavior is exemplified in Figure [TJ where panels (a) and (c) show the time evolution of 
the GALI2 index for two FPU trajectories started in the vicinity of 2-tori of the = 32 
and A^ = 128 systems, when /3 = 0.1 and tmax = 10''. In both cases, the energy satisfies 
E > Ec, as the exponential fall-off of the GALI2 index is already observed at t = 10''. 
However, a simple inspection of Figs. 03, d clearly reveals that the exponential localization 
of the energy persists in the Fourier space of both systems. In fact we have found that the 
exponential localization persists for energies much larger than Ec, but for timescales which 
become smaller as E increases. 



B. A heuristic argument on exponential stability 



The results of the previous subsection confirm the observation made on the basis of numer- 
ical experiments in lo|, that a result establishing exponential stability for FPU trajectories 
close to g-tori should be possible. This has already been discussed in the introduction, where 
we mentioned a partial result in this direction obtained recently in [l^]. It was pointed out, 
however, that so far the main obstruction to precise analytical statements lies in the bad 
dependence of all existing estimates on A^. 

Thus, in closing our paper, we would like to offer a heuristic argument on how the 
property of exponential localization of the energy for particular FPU trajectories could serve 
as a basis for further improvement of rigorous results. For definiteness, we refer below to 
the case of FPU-trajectories for which the question of exponential stability can be examined 
in the framework of a variant of Birkhoff's method due to Giorgilli [2^, based on the direct 



calculation of approximate integrals of motion without the use of normal forms [40|, |41|, 142 1. 

The trajectories we are referring to have initial conditions close to a particular family of 
g-breather solutions considered in 38|, 1391], whose 'seed mode' go varies proportionally to A^. 
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We thus take the seed mode to vary as qo = {N + l)/2/i, where fi can be any large number 
that is a power of two, while is a power of two minus 1, so that no commensurabilities 
exist in the set of unperturbed FPU frequencies Qg, q = 1,...N. The exact g-breather 
solution oscillates with just one frequency and thus lies on a one-dimensional torus. However, 
perturbing this trajectory, we obtain solutions lying on /i— dimensional tori involving the 
modes q = {2j — l)qo, j = 1,2, Furthermore, the energy localization profile for the 

g-breather solution is the same as that of a s-dimensional g-torus solution with s = qo given 
by Eq.([35]), where M = go/^V - 1/2/i. 

Let us assume that this localization profile holds also initially for the FPU-trajectory 
that is a perturbation of the exact g-breather solution. Our aim then is to demonstrate that 
when one establishes estimates of exponential time of stability T ~ exp(l/£:'^) for the energy 
localization profile of this type of FPU-trajectories, the exponent c does not depend on the 
number of degrees of freedom (i.e. is independent of /i), despite the diophantine character 
of the unperturbed frequencies Qq of the fj, modes g = (2j — l)go, j = 1, 2, 

To this end, we recall that, according to 2^, the question of the long-term stability of 
particular trajectories can be examined by constructing approximate integrals of motion in 
the form of truncated series starting with harmonic energies 

^^ = ^f + + (37) 

where is an rth degree polynomial in the canonical variables (gi, . . . , gAr,pi, . . . ,pjv) 
defined by 



and (^'^p = Ej = i^ljPjqj. The functions are determined recursively, by solving formally 
the Poisson bracket condition + + ■■■1H2 + H4} = 0, demanding that the series 
$j be an integral. Here H2 and H4 are the quadratic and quartic terms of the /3-FPU 
Hamiltonian written in the variables {qj,pj). The order by order solution is then found via 
the so-called homological equation 

H,} + {<l>f , H2} = 0, r = 4, 6, ... (39) 

which has the same algebraic structure as the homological equation of Birkhoff 's normaliza- 
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tion scheme. The solution reads 

(r) _ _ V- hm,nq P . . 

^ ^ (m-n)-Q ^ ^ 

m+n — r 

where use is made of the compact notation q™p" = q^^q^^ ■ ■ ■ q^^'pTT>2' ■ ■ -Pn^ ^ ™ = 
{■1711,1712, . . . ,171^), m = mi + m2 + . . . + 171^, (similarly for n), and hm,n are the polynomial 
coefficients of the Poisson bracket {$^'^ ^^7^4}. 

It is well known that the series ( !37l) is asymptotic. We thus consider the rth order finite 
truncation $,> = + (^^f' + . . . + (^^p representing an approximate integral of motion 
whose time variation is given by 

^ = r,,^{^^;\h,} . (41) 

The quantity Rj^r is called the remainder function. The asymptotic character of the series 
implies that the size of Rj^r decreases initially as r increases, up to an optimal order r = Vopt 
beyond which the size of the remainder increases with r, becoming ultimately divergent as 
r ^ 00. Recursive application of the homological equation fl39l) (see [20|, |2J] for details) 
yields that the size of the remainder at order r can be estimated by 

ir-2)\\R,,.-2\\ (r-2)(r-4)||i?,,,_4|| jr - 2)\\\\R,,,\\ 

1 1 -thj^r 1 1 ~ ~ ~ . . . ~ [^■^J 

where | |-Rj,r| | denotes the absolute sum of the polynomial coefficients n for all (m, n) with 
r and Ofc, k = 4, 6, r — 2, denotes the minimum of all divisors (m — n) • fi, 
= k at order k. The pre-factors (r — 2), (r — 4)... come from the derivatives 
of the functions ^^J^ ^^J^ ^\---, with respect to the canonical variables (q, p), appearing 
in Poisson brackets due to recursive application of the homological equation, which yield 
factors |m|, |n|, both of order 0{r — 2), 0{r — 4),... within the functions ^^^~'^\ $ 



m - 


- n 


m - 


F |n| 



(r-2) ^(r-4) 

, . . . , 



(4) 



Our heuristic argument now goes as follows: Since in the considered trajectories the 
total number of participating modes is equal to /i, and is a power of two minus one, the 
associated divisors satisfy a diophantine condition of the form 

|1 • f^l > ^ for all 1 e Z^, |I| ^ (43) 
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R,,r\\-{r-2)\\^^'C/[l] (44) 



with r > /i — 1. The estimate r ~ /i holds for large. One then has 
^%^(r — 2)!!^'', which, upon substitution in (H2l) leads to: 

r/2 

where the estimate 1 1-^41 1 = 0{j3/N) is taken into account and C is a 0(1) positive con- 
stant. These estimates, based on small divisors, are standard and lead by themselves to 
no improvement as far as the dependence of the asymptotic character on fi is concerned. 
However, an improvement can be achieved if we also take into account the numerators of 
the remainder series, which for FPU-trajectories possess the exponential localization profile 
shown in Eq.( lMl) . In other words, the size of the remainder depends also on the size of the 
monomials 

where C,k stands for either q^. or p^, since each term of the remainder consists of one such 
monomial multiplied by a coefficient bounded by an estimate of the form (jH]). In view of 
Eq.( !38l) . one has ~ eI^'^. Thus, taking into account the form of the profile (l34|) . the 

size of the above monomial at the rth order of normalization can be estimated as: 

|. I £2 .„xf±i s-i^+S2 + ...+st, 2 \S2+2s3 + ... + {/i-l)s^ 

|A^| ~ (^l) 2 {E2) 2 ...{E^) 2 E 2 [PfJ' £) 

where Si -|- S2 + • • . + = r. For the leading terms (with smallest divisors) in the series the 
exponents Si, S2, typically take values such that the estimate S2 + 2S3 + . . . + (/i — l)s^ ~ 
fir/ 2 holds. Thus 

|A^| ~ VJE'l^e^'l^ (45) 

where is a new constant. Combining now the estimates (145|) and fH4|) yields a final 
estimate for the size of the remainder of the form 

~ (r - 2)\\^^^CJVJ ' ~ ~ B,{{r - 2y)^e^^+^>/^ , (46) 

The optimal order of truncation is found by taking the logarithm of ( 146|) . using Stirling's 
formula logra! ^ nlogra — n, as well as r — 2 ^ r (for r large), using rflog ~ 0. We 

thus find ropt ~ 1/e, whence the optimal value of the remainder is 

ll^op*||~expf-^) (47) 



2e 
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i.e. the time variations of the integrals ^j,ropt exponentially small in 1/e. This yields an 
estimate of exponential stability of the form T ~ exp(l/£:^), where c = 1, i.e. the number 
of degrees of freedom no longer appears in the exponent c. Compared to general estimates 
yielding c ~ it is seen that the removal of the dependence of c on /i was possible 

thanks to the assumed exponential scaling of the energy localization profile, that allows 
for the estimate (H5l) . In turn, since this localization profile is the same as for g-tori, the 
above analysis is suggestive of the usefulness of exponential localization in g-space in order 
to obtain improved estimates in the case of g-tori solutions as well. However, in the lack 
of a rigorous demonstration, it remains an open question whether such strategy can lead to 
estimates of a stretched exponential law for the dependence of ||-Ropt|| on e, as evidenced 
also to a limited extent by the numerical experiments of [lo| . 



IV. CONCLUSIONS 



The main conclusions of the present study can be summarized as follows: 

1) We introduced the concept of g-tori in FPU lattices, which represent a generalization 
of the concept of g-breathers. The g-tori have low dimensionality s << N, and arise from 
the continuation of motions with s independent frequencies of the unperturbed problem. 

2) We explicitly calculated FPU solutions lying on g-tori by employing the method of 
Poincare-Lindstedt series. Based on estimates of the leading terms of these series, we pro- 
vided a theoretical law yielding the average exponential localization of the energy in Fourier 
space for solutions on g-tori. Furthermore, a proposition was proved which explains how 
different groups of modes are excited at consecutive orders of the perturbation theory. The 
most important conclusion from this analysis is that, if the fraction s/N is kept constant, it 
appears that the localization profile depends on the specific energy e = E/N of the system, 
i.e., is independent of A^. Some numerical evidence is provided in support of the above 
theoretical analysis, but further numerical work is necessary, in order to clarify the extent 
of its validity. 

3) We explored numerically the relevance of g-tori of dimension s to the dynamics of 
FPU-trajectories started nearby, by exciting s modes only. The localization laws found 
analytically for g-tori accurately describe the localization of energy in Fourier space for 
the FPU-trajectories as well. We also gave numerical evidence of the existence of two 
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regimes separated by a critical energy value. Below this energy, the localization profiles 
E[q) oc exp(— 6g) have a slope b depending logarithmically on e, while the fraction M = s/N 
of modes sharing the energy is constant. Furthermore, beyond this critical energy the slope 
of the localization profile tends to stabilize, and M tends asymptotically to the law M oc e^^'^, 
as g-tori of progressively higher dimension begin to describe the dynamics of the numerical 
FPU-trajectories. 

4) We examined the stability of g-tori using a numerical criterion provided by the GALI 
indices j^^l , and provided numerical evidence demonstrating that the localization in Fourier 
space persists for energies well above the threshold value at which the underlying g-tori turn 
from stable to unstable. 

5) Finally, we provided a heuristic argument suggesting that the exponential energy 
localization in g-mode space implies a particular analytical structure of the Birkhoff series 
obtained for associated FPU-trajectories, that could lead to improved estimates of the long 
term stability in the spirit of Nekhoroshev theory. In the example of FPU trajectories 
started close to g-breather solutions exhibiting the same energy localisation as for g-tori, 
our arguments suggest that in stability estimates of the form T ~ exp(l/e^) it is possible 
to remove the bad dependence of c on the number of degrees of freedom of the problem. 
Of course, further study, substantiated by numerical experiments, is needed before rigorous 
statements are available on this issue. 
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APPENDIX A: SEQUENCE OF MODE EXCITATIONS 

For all g 7^ gj, i = 1, . . . s, the equations yielding the Poincare-Lindstedt series terms Q^g'' 
at first order are 

N 

Q^^ + ^lQ^i^ = - Yl ^,^i^ACgi„,nQTQ^^^Q'^^^ . (Al) 

i,m,n=l 

The r.h.s. of the above equations is different from zero if Cgimn 7^ 0. In view of the definition 
of Cgimn, non-zero solutions of Eqs. flAip containing no new frequencies fig are possible 
for the values of q satisfying either q = q^^^ =\ ±/±m±n |= mi with l,m,n E {qi, q2, qs}, 
when 1 < mi < A^, or g^^^ =| 2{N + 1)- \ ±1 ±m ±n \ \, when + 2 < mi < 2A^ + 1 or 
2A^+3 < mi < 3A^. Both cases are given by Eq. fl26|) . for A = and A = 1 respectively. Thus, 
the proposition holds for A; = 1. Assuming it to be true at order k — 1, and using Eq. (l3T!) . 
one finds for the kth order that solutions introducing no new frequencies are possible for the 
modes satisfying either q = q^''^ =\ ±/±m±n |, 1 < q^'^^ < N,oi q = q^^^ =| 2(A^ + 1)— mi |, 
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1 < g(^) < A^, where 

/ = |2A„^(A^+ 1) -m„J, m„^ = ± ± ... ± g^j^^^J, 
m = \2Xn^{N + 1) -rrin^l, m^^ = jg^-, ± g^-^ ± ... ± g,. |, A„2 



n = |2A„3(A^ + 1) -m 



na I ; 



.2(iV + l). 

' rrin,, + N 
.2(iV + l) 

m„3 + iV 



(A2) 



.2(Ar + l) 

with ji, j2ni+i,?^i, •••,?^2n2+i,^i, •••,^2n3+i ^ {l,...,s}, ni + n2 + ^3 = - 1 • Taking the last 
relation as well as all possible sign combinations in the sum ±/ ± m ± n into account, the 
permitted modes at fcth order are given by equations of the form 



|2(±A„, ± A„, ± A„3 ±g){N + l) + (±g,, ± g,, ± . . . ± g,,,^^; 



(A3) 



provided that 1 < g'^'^^ < A^, where = or 1, and ii,i2, ■ ■ ■ ,i2k+i G {1, . . . , s}. After a 
possible sign reversal within | • |, not affecting the absolute value, the expression 



|2(±A„, ± Xn, ± A„3 ±g){N+l) + (ig^, ± g^, ± . . . ± q^^^^^] 



always resumes the form 



g^'^) = |2A(A^+ 1) -mfcl 



(A4) 



where ruk = |gii ± gij ± • • • ± giafe+J and A is an integer number. However, by the second 
restriction of (lA3p . namely 1 < g*^'^^ < A^, one necessarily has that A = [(m^ + A^)/2(A^ + 1)], 
which concludes the proof of the proposition. 



APPENDIX B: ESTIMATES ON LOCALIZATION PROFILES 

The products <5^"^])Q^"^2)<5^"4\ Eg- (131]) give rise to trigonometric terms of the form 



Jni) ^(n2) ^(na) 
'g{ni)Vg{n2)Vg{' 

tti G Z, \ai\ 



Q2^l)Qj^(ll)Qg(ll) cos(aiu;i + a2U2 + ... + asU},)t : 

lazl + ... + laJ <2k + l 



(Bl) 



on the r.h.s. of Eq. (l3T|) . For each trigonometric term of the form (IBip . and A; > 1, Eq.(l3T 



introduces a divisor to the solution for Q\l), namely 



(fe) 



(B2) 
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In view of Eg. (1220 . the smallest divisors are those satisfying 



For such divisors, one has 



while, for the sum X]y=i '^jJ^ has the inequality 



O {{nsk/Nf) 



j=i j=i j=i j=i 



in view of which the estimate 



j=i 



4 I 



(B3) 



holds. On the other hand, using condition (IBSp as well as Eq. (!28|) one obtains 

s s ^2 ^ 

j=i j=i j=i 
Combining the last two expressions we find 

i=i 

in view of which (combined with fl28l) ) one finally arrives at an estimate for the size of 
divisors 



UJ (k) 



-E 



7r3s3(2A; + l)[(2A;+ 1)2-1] 27rs(2fc + 1) 



24(iV+ 1)3 



(iV + 1) 



(B4) 



12(A^ + 1)4 

Returning to Eq. (l3T!) . for fixed values of ni, n2 and ns, there are at most s triplets 
g{"2)^ g(n3)) satlsfylug g("2)^ g{n3)) g . Usiug this fact as well as the esti- 

mates (lB4p and (1281) . one obtains for the norms of the various terms the estimate: 

12(iV + l)4 7rs(2A; + l) 



■X 
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fc-1 



^ (2ni + l)(2n2 + l)(2n3 + l)||Q(;\) II HQ^^^H \\Q^2l 



"1,2,3=0 

ri]^ +712 +^13 — fc — 1 

or 



\\Q%\\-J^^TW ^ (2m + l)(2n, + l)(2n3 + l)||gj;\\|| Hg^-^H ||Q^:^3)II(B5) 

Denoting by A*^^) the average size of the oscillations of all the modes Eq. flBSp takes the 
form 

12s '"'^ 

-— — ^ V (2ni + l) (2^2 + 1) (2^3 + 1)A("^)A("^M("«) . (B6) 



By induction it now follows that 

\n /\2n+l 



(B7) 

2n + l ^ ^ 

Indeed, assuming fiBTl) to be true for the amplitudes A^''^\ A^^^), ^("3) Eq.(lB6]), it follows 
that 

^ (3s/2)^-H2.Ar^ ^ . ^ 

(2A;+1)3 2^ 

"1,2,3=0 

n-j^ +ri2 +^1-3 — fc — 1 

{3s/2)^-^12sAl''+^ k{k + l) {3s/2)''-'l2sAl^+^ (2A; + l)^ {3s/2)''Al^+^ 



(2A; + 1)3 2 (2A; + 1)3 8 2k + 1 

which demonstrates the validity of Eq.( l32|) in the text, with the constant C having the 
specific value C = 3/2. 
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